#include "Convoluter.h"

Convoluter::Convoluter(ImportData* d)
{	
	this->data = d;
}

void Convoluter::initialize2d(){
	//Set Filter Size
	const int size = 3;
	// Initialize the structures
	dxFilter.resize(size,size);
	dyFilter.resize(size,size);
	dxdxFilter.resize(size,size);
	dxdyFilter.resize(size,size);
	dydyFilter.resize(size,size);

	/// Initialize the filters //original 2d
    dxFilter << -1,0,1, -2,0,2, -1,0,1;
    dyFilter << 1,2,1, 0,0,0, -1,-2,-1;
    dxdxFilter << 1,-2,1,2,-4,2,1,-2,1;
    dxdyFilter << 1,0,-1,0,0,0,-1,0,1;
    dydyFilter << 1,2,1,-2,-4,-2,1,2,1;

	/// Initialize the filters //Match3d
    //dxFilter << -2,0,2, -4,0,4, -2,0,2;
    //dyFilter << 2,4,2, 0,0,0, -2,-4,-2;
    //dxdxFilter << 2,-4,2,4,-8,4,2,-4,2;
    //dxdyFilter << 2,0,-2,0,0,0,-2,0,2;
    //dydyFilter << 2,4,2,-4,-8,-4,2,4,2;

	/*cout << "CONVOLUTER FILTERS: " << endl;
	cout << dxFilter << '\n' << endl;
	cout << dyFilter << '\n' << endl;
	cout << dxdxFilter << '\n' << endl;
	cout << dxdyFilter << '\n' << endl;
	cout << dydyFilter << '\n' << endl;*/
}

void Convoluter::initialize3d(){
	//Set Filter Size
	const int size = 3;
	// Initialize the structures
	dxFilter3d.resize(size); //depth of filter
	dyFilter3d.resize(size); //depth of filter
	dzFilter3d.resize(size); //depth of filter

	dxdxFilter3d.resize(size); //depth of filter
	dydyFilter3d.resize(size); //depth of filter
	dzdzFilter3d.resize(size); //depth of filter
	dxdzFilter3d.resize(size); //depth of filter
	dxdyFilter3d.resize(size); //depth of filter
	dydzFilter3d.resize(size); //depth of filter

	for(int i = 0;i<size;++i){
		dxFilter3d[i].resize(size,size);
		dyFilter3d[i].resize(size,size);
		dzFilter3d[i].resize(size,size);

		dxdxFilter3d[i].resize(size,size);
		dydyFilter3d[i].resize(size,size);
		dzdzFilter3d[i].resize(size,size);
		dxdzFilter3d[i].resize(size,size);
		dxdyFilter3d[i].resize(size,size);
		dydzFilter3d[i].resize(size,size);
	}
		
	/// Initialize the filters
	dxFilter3d[0] << -1,0,1,-2,0,2,-1,0,1; //bottom
	dxFilter3d[1] << -2,0,2,-4,0,4,-2,0,2; //middle
	dxFilter3d[2] << -1,0,1,-2,0,2,-1,0,1; //top

	dzFilter3d[0] << 1,2,1,0,0,0,-1,-2,-1; //bottom
	dzFilter3d[1] << 2,4,2,0,0,0,-2,-4,-2; //middle
	dzFilter3d[2] << 1,2,1,0,0,0,-1,-2,-1; //top

	dyFilter3d[0] << -1,-2,-1,-2,-4,-2,-1,-2,-1; //bottom
	dyFilter3d[1] << 0,0,0,0,0,0,0,0,0; //middle
	dyFilter3d[2] << 1,2,1,2,4,2,1,2,1; //top

	//dxdxFilter3d,dydyFilter3d,dzdzFilter3d,dxdzFilter3d,dxdyFilter3d,dydzFilter3d;

	dxdxFilter3d[0] << 1,-2,1,2,-4,2,1,-2,1; //bottom
	dxdxFilter3d[1] << 2,-4,2,4,-8,4,2,-4,2; //middle
	dxdxFilter3d[2] << 1,-2,1,2,-4,2,1,-2,1; //top

	dydyFilter3d[0] << 1,2,1,2,4,2,1,2,1; //bottom
	dydyFilter3d[1] << -2,-4,-2,-4,-8,-4,-2,-4,-2; //middle
	dydyFilter3d[2] << 1,2,1,2,4,2,1,2,1; //top

	dzdzFilter3d[0] << 1,2,1,-2,-4,-2,1,2,1; //bottom
	dzdzFilter3d[1] << 2,4,2,-4,-8,-4,2,4,2; //middle
	dzdzFilter3d[2] << 1,2,1,-2,-4,-2,1,2,1; //top

	//dxdyFilter3d[0] << 2,0,-2,4,0,-4,2,0,-2; //bottom
	//dxdyFilter3d[1] << 0,0,0,0,0,0,0,0,0; //middle
	//dxdyFilter3d[2] << -2,0,2,-4,0,4,-2,0,-2; //top

	//dxdzFilter3d[0] << 2,0,-2,4,0,-4,2,0,-2; //bottom
	//dxdzFilter3d[0] << 0,0,0,0,0,0,0,0,0; //middle
	//dxdzFilter3d[0] << -2,0,2,-4,0,4,-2,0,-2; //top

	//dydzFilter3d[0] << 2,0,-2,4,0,-4,2,0,-2; //bottom
	//dydzFilter3d[0] << 0,0,0,0,0,0,0,0,0; //middle
	//dydzFilter3d[0] << -2,0,2,-4,0,4,-2,0,-2; //top
}

void Convoluter::computeDerivatives(int slice, MatrixXd& outDx,MatrixXd& outDy,
											MatrixXd& outDxDx,MatrixXd& outDxDy,MatrixXd& outDyDy) const

{	
	cout << "CONVOLUTING DATA OVER SLICE " << slice << " THAT HAS WxH: " << this->data->getWidth() << "," << this->data->getHeight() << "\n" << endl;
	//Define Indata
	DataPoint** in;
	in = new DataPoint*[this->data->getWidth()];
	//Initialize Indata
	for(int i = 0; i < this->data->getWidth(); ++i){
		in[i] = new DataPoint[this->data->getHeight()];
	}

	//Designate the 2D slice of the DataVolume That is Convoluted
	int x, y;
    for(x = 0; x < this->data->getWidth(); ++x){
        for(y = 0; y < this->data->getHeight(); ++y){
			in[x][y] = this->data->getVolume()[slice][y][x];
		}
	}

    // Initialize all the 5 output variables so that they are of the same size as in
    outDx.resize(this->data->getWidth(),this->data->getHeight());
	outDy.resize(this->data->getWidth(),this->data->getHeight());
	outDxDx.resize(this->data->getWidth(),this->data->getHeight());
	outDxDy.resize(this->data->getWidth(),this->data->getHeight());
	outDyDy.resize(this->data->getWidth(),this->data->getHeight());
    
    // compute the derivatives 
    convolute(in, this->dxFilter, outDx);
    convolute(in, this->dyFilter, outDy);
    convolute(in, this->dxdxFilter, outDxDx);
    convolute(in, this->dxdyFilter, outDxDy);
    convolute(in, this->dydyFilter, outDyDy);

	delete in; //Remove excess memory used by in

} /* computeDerivatives() */

void Convoluter::computeDerivatives3d(std::vector<MatrixXd>& outDx,std::vector<MatrixXd>& outDy,std::vector<MatrixXd>& outDz,std::vector<MatrixXd>& outDxDx,
									  std::vector<MatrixXd>& outDyDy,std::vector<MatrixXd>& outDzDz,std::vector<MatrixXd>& outDxDy,
									  std::vector<MatrixXd>& outDxDz,std::vector<MatrixXd>& outDyDz) const{
	// Initialize the output variable so that it is of the same size as volume
	outDx.resize(this->data->getDepth()); //dx
	outDy.resize(this->data->getDepth()); //dy
	outDz.resize(this->data->getDepth()); //dz
	
	outDxDx.resize(this->data->getDepth()); //dxdx
	outDyDy.resize(this->data->getDepth()); //dydy
	outDzDz.resize(this->data->getDepth()); //dzdz
	outDxDy.resize(this->data->getDepth()); //dxdy / dydx
	outDxDz.resize(this->data->getDepth()); //dxdz / dzdx
	outDyDz.resize(this->data->getDepth()); //dydz / dzdy

	for(int i = 0;i<outDz.size(); ++i){
		outDz[i].resize(this->data->getWidth(),this->data->getHeight()); //x,y size
		outDx[i].resize(this->data->getWidth(),this->data->getHeight()); //x,y size
		outDy[i].resize(this->data->getWidth(),this->data->getHeight()); //x,y size

		outDxDx[i].resize(this->data->getWidth(),this->data->getHeight()); //x,y size
		outDyDy[i].resize(this->data->getWidth(),this->data->getHeight()); //x,y size
		outDzDz[i].resize(this->data->getWidth(),this->data->getHeight()); //x,y size
		outDxDy[i].resize(this->data->getWidth(),this->data->getHeight()); //x,y size
		outDxDz[i].resize(this->data->getWidth(),this->data->getHeight()); //x,y size
		outDyDz[i].resize(this->data->getWidth(),this->data->getHeight()); //x,y size
	}
	// compute the derivatives 
	convolute3D(this->dxFilter3d,outDx);
	convolute3D(this->dyFilter3d,outDy);
    convolute3D(this->dzFilter3d,outDz);

	convolute3D(this->dxdxFilter3d,outDxDx);
	convolute3D(this->dydyFilter3d,outDyDy);
	convolute3D(this->dzdzFilter3d,outDzDz);

	////dxdy
	//convolute3D(this->dxFilter3d,outDxDy);
	//convolute_second_3D(outDxDy,this->dyFilter3d,outDxDy);
	convolute_second_3D(outDx,this->dyFilter3d,outDxDy);
	////dxdz
	//convolute3D(this->dxFilter3d,outDxDz);
	convolute_second_3D(outDx,this->dzFilter3d,outDxDz);
	////dydz
	//convolute3D(this->dyFilter3d,outDyDz);
	convolute_second_3D(outDy,this->dzFilter3d,outDyDz);
}

void Convoluter::convolute3D(const std::vector<MatrixXd>& filter,std::vector<MatrixXd>& out) const {
	//in: this->data->getVolume();
	const int F_W_MIN = -((filter[0].cols() - 1) / 2);
    const int F_W_MAX = filter[0].cols() / 2;
    const int F_H_MIN = -((filter[0].rows() - 1) / 2);
    const int F_H_MAX = filter[0].rows() / 2; 
	const int F_D_MIN = -((filter.size() - 1) / 2);
	const int F_D_MAX = filter.size() / 2;

	/*cout << "F_W: " << F_W_MIN << "," << F_W_MAX << endl;
	cout << "F_H: " << F_H_MIN << "," << F_H_MAX << endl;
	cout << "F_D: " << F_D_MIN << "," << F_D_MAX << endl;*/

	int x,y,z,i,j,k;
	 for(x = 0; x < this->data->getWidth(); ++x){
        for(y = 0; y < this->data->getHeight(); ++y){
			for(z = 0; z < this->data->getDepth(); ++z){
				if(!this->data->getVolume()[z][y][x].isEmpty()){ // if the point is valid/not empty
				 out[z](x,y) = 0;
				 // and do the convolution at (x, y) for each z -1,0,1
					for(i = F_W_MIN; i <= F_W_MAX; ++i){ //x
						for(j = F_H_MIN; j <= F_H_MAX; ++j){ //y 
							for(k = F_D_MIN; k <= F_D_MAX; ++k){ //z
								out[z](x,y) += filter[k-F_D_MIN](i - F_W_MIN, j - F_H_MIN) * getValueAt3d(x + i, y + j,z + k);
							}//k
						}//j
					}//i
				}//if !empty
			} //z
		} //y
	 } //x

}

void Convoluter::convolute_second_3D(std::vector<MatrixXd>& in,const std::vector<MatrixXd>& filter,std::vector<MatrixXd>& out) const {
	//in: this->data->getVolume();
	const int F_W_MIN = -((filter[0].cols() - 1) / 2);
    const int F_W_MAX = filter[0].cols() / 2;
    const int F_H_MIN = -((filter[0].rows() - 1) / 2);
    const int F_H_MAX = filter[0].rows() / 2; 
	const int F_D_MIN = -((filter.size() - 1) / 2);
	const int F_D_MAX = filter.size() / 2;

	/*cout << "F_W: " << F_W_MIN << "," << F_W_MAX << endl;
	cout << "F_H: " << F_H_MIN << "," << F_H_MAX << endl;
	cout << "F_D: " << F_D_MIN << "," << F_D_MAX << endl;*/

	int x,y,z,i,j,k;
	 for(x = 0; x < this->data->getWidth(); ++x){
        for(y = 0; y < this->data->getHeight(); ++y){
			for(z = 0; z < this->data->getDepth(); ++z){
				if(!this->data->getVolume()[z][y][x].isEmpty()){ // if the point is valid/not empty
				 out[z](x,y) = 0;
				 // and do the convolution at (x, y) for each z -1,0,1
					for(i = F_W_MIN; i <= F_W_MAX; ++i){ //x
						for(j = F_H_MIN; j <= F_H_MAX; ++j){ //y 
							for(k = F_D_MIN; k <= F_D_MAX; ++k){ //z
								out[z](x,y) += filter[k-F_D_MIN](i - F_W_MIN, j - F_H_MIN) * getValueAt_second_3d(in,x + i, y + j,z + k);
							}//k
						}//j
					}//i
				}//if !empty
			} //z
		} //y
	 } //x

}

double Convoluter::getValueAt_second_3d(std::vector<MatrixXd>& in,int x, int y,int z) const
{
    // Here you check what value you should have at a certain location:
    // If (x, y) is valid and in the domain, you return in[x][y]
    // If (x, y) is not valid or outside of the domain (lesser than 0 or 
    // greater than width or height) you do the mirroring
	//cout << "Asking For Value at: " << x << "," << y << endl;
	
	if(z < 0){
		z = -z;
	}
	else if(z > (this->data->getDepth()-1)){ // mirror z if necessary
		z = (this->data->getDepth()-1) - 1 - (z - this->data->getDepth());
	}
	if(x < 0){ // mirror x if necessary
		//cout << "mirror x as negative" << endl;
        x = -x;
	}
	else if(x > (this->data->getWidth()-1)){
		//cout << "mirror x as out of bounds" << endl;
        x = (this->data->getWidth()-1) - 1 - (x - this->data->getWidth());
        // or x = 2 * in.getWidth() - x - 1
	}
	if(y < 0){ // mirror y if necessary
		//cout << "mirror y as negative" << endl;
        y = -y;
		//cout << "mirror y as negative done" << endl;
	}
	else if(y > (this->data->getHeight()-1)){
		//cout << "mirror y as out of bounds" << endl;
        y = (this->data->getHeight()-1) - 1 - (y - this->data->getHeight());
        // or y = 2 * in.getHeight() - y - 1
	}
	///cout << "getting value for point at " << x << "," << y << endl;
	if(this->data->getVolume()[z][y][x].isEmpty()){ // if the point is not valid...
        // mirroring in this case is a bit too complex, unless we see that 
        // it does not work, we can discard this point
		//cout << "getvalue finished == 0" << endl;
		//cout << "point not valid" << endl;
        return 0;
    }
    // If the point is valid, we return the value at that position
	//cout << "returning: " << in[z](x,y) << endl;;
	return in[z](x,y);

} /* getValueAt() */

double Convoluter::getValueAt3d(int x, int y,int z) const
{
    // Here you check what value you should have at a certain location:
    // If (x, y) is valid and in the domain, you return in[x][y]
    // If (x, y) is not valid or outside of the domain (lesser than 0 or 
    // greater than width or height) you do the mirroring
	//cout << "Asking For Value at: " << x << "," << y << endl;
	
	if(z < 0){
		z = -z;
	}
	else if(z > (this->data->getDepth()-1)){ // mirror z if necessary
		z = (this->data->getDepth()-1) - 1 - (z - this->data->getDepth());
	}
	if(x < 0){ // mirror x if necessary
		//cout << "mirror x as negative" << endl;
        x = -x;
	}
	else if(x > (this->data->getWidth()-1)){
		//cout << "mirror x as out of bounds" << endl;
        x = (this->data->getWidth()-1) - 1 - (x - this->data->getWidth());
        // or x = 2 * in.getWidth() - x - 1
	}
	if(y < 0){ // mirror y if necessary
		//cout << "mirror y as negative" << endl;
        y = -y;
		//cout << "mirror y as negative done" << endl;
	}
	else if(y > (this->data->getHeight()-1)){
		//cout << "mirror y as out of bounds" << endl;
        y = (this->data->getHeight()-1) - 1 - (y - this->data->getHeight());
        // or y = 2 * in.getHeight() - y - 1
	}
	///cout << "getting value for point at " << x << "," << y << endl;
	if(this->data->getVolume()[z][y][x].isEmpty()){ // if the point is not valid...
        // mirroring in this case is a bit too complex, unless we see that 
        // it does not work, we can discard this point
		//cout << "getvalue finished == 0" << endl;
		//cout << "point not valid" << endl;
        return 0;
    }
    // If the point is valid, we return the value at that position
	return this->data->getVolume()[z][y][x].getFTLE();

} /* getValueAt() */


//template <class Data, FilterT>
void Convoluter::convolute(DataPoint** in, const MatrixXd& filter,
        MatrixXd& out) const 
{
    const int F_W_MIN = -((filter.cols() - 1) / 2);
    const int F_W_MAX = filter.cols() / 2;
    const int F_H_MIN = -((filter.rows() - 1) / 2);
    const int F_H_MAX = filter.rows() / 2; 

    // If the filter is 3x3 the previous values are -1, 1, -1, 1
    // If the filter is 4x4 the previous values are -1, 2, -1, 2
    // If the filter is 5x5 the previous values are -2, 2, -2, 2
    // and so on

    int x, y, i, j;
    for(x = 0; x < this->data->getWidth(); ++x){
        for(y = 0; y < this->data->getHeight(); ++y){
			//cout << "first for: " << x << "," << y << endl;
            //if(in.isValid(x, y){ // if the point is valid
			if(!in[x][y].isEmpty()){ // if the point is valid/not empty
                out(x,y) = 0; // init the output variable
                // and do the convolution at (x, y)
                for(i = F_W_MIN; i <= F_W_MAX; ++i){
                    for(j = F_H_MIN; j <= F_H_MAX; ++j){  
						//cout << "getting value at " << (x+1) << "," << (y+j) << endl;
                        out(x,y) += filter(i - F_W_MIN, j - F_H_MIN) * getValueAt(in, x + i, y + j);
                    } /* for j */
                } /* for i */
				//cout << "outvalue at " << x << "," << y << " = " << out(x,y) << endl;
            } /* if (x, y) is valid */
        } /* for y */
    } /* for x */
} /* convolute() */



double Convoluter::getValueAt(DataPoint** in, int x, int y) const
{
    // Here you check what value you should have at a certain location:
    // If (x, y) is valid and in the domain, you return in[x][y]
    // If (x, y) is not valid or outside of the domain (lesser than 0 or 
    // greater than width or height) you do the mirroring
	//cout << "Asking For Value at: " << x << "," << y << endl;

	if(x < 0){ // mirror x if necessary
		//cout << "mirror x as negative" << endl;
        x = -x;
	}
	else if(x > (this->data->getWidth()-1)){
		//cout << "mirror x as out of bounds" << endl;
        x = (this->data->getWidth()-1) - 1 - (x - this->data->getWidth());
        // or x = 2 * in.getWidth() - x - 1
	}
	if(y < 0){ // mirror y if necessary
		//cout << "mirror y as negative" << endl;
        y = -y;
		//cout << "mirror y as negative done" << endl;
	}
	else if(y > (this->data->getHeight()-1)){
		//cout << "mirror y as out of bounds" << endl;
        y = (this->data->getHeight()-1) - 1 - (y - this->data->getHeight());
        // or y = 2 * in.getHeight() - y - 1
	}
	///cout << "getting value for point at " << x << "," << y << endl;
	if(in[x][y].isEmpty()){ // if the point is not valid...
        // mirroring in this case is a bit too complex, unless we see that 
        // it does not work, we can discard this point
		//cout << "getvalue finished == 0" << endl;
		//cout << "point not valid" << endl;
        return 0;
    }
    // If the point is valid, we return the value at that position
	return in[x][y].getFTLE();

} /* getValueAt() */


Convoluter::~Convoluter(void)
{
}
